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Abstract 

The classical theory of ion beam sputtering predicts the instability of a flat surface to uniform 
ion irradiation at any incidence angle. We relax the assumption of the classical theory that the 
average surface erosion rate is determined by a Gaussian response function representing the effect 
of the collision cascade and consider the surface dynamics for other physically-motivated response 
functions. We show that although instability of flat surfaces at any beam angle results from 
all Gaussian and a wide class of non-Gaussian erosive response functions, there exist classes of 
modifications to the response that can have a dramatic effect. In contrast to the classical theory, 
these types of response render the flat surface linearly stable, while imperceptibly modifying the 
predicted sputter yield vs. incidence angle. We discuss the possibility that such corrections underlie 
recent reports of a "window of stability" of ion-bombarded surfaces at a range of beam angles for 
certain ion and surface types, and describe some characteristic aspects of pattern evolution near the 
transition from unstable to stable dynamics. We point out that careful analysis of the transition 
regime may provide valuable tests for the consistency of any theory of pattern formation on ion 
sputtered surfaces. 

PACS numbers: 68.49. Sf, 81.65.Cf, 81.16.Rf 
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I. INTRODUCTION 



Uniform ion beam sputter erosion of a solid surface often causes a spontaneously- arisin, 



topographic pattern in the surface t opogra phy 
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30, l3l|, I33, l33| , that can take the form of 



a one-dimensional corrugation or a two-dimensional array of dots with typical length scales 

of 10 nm. Periodic self-organized patterns with wavelength as small as 15 nm 

have stimulated interest in this method as a means of nanofabrication at sub-lithographic 



length scales 



26|. Because the characteristic scale of the patterns can be three orders of 



magnitude larger than the characteristic penetration depth of ions into a solid surface, the 
patterns result from a nontrivial interplay between the sputter erosion on one hand and 
surface relaxation mechanisms on the other hand. 

The present understanding of sputter morphology evolution originates in the Sigmund 
theory of sputtering [s^. Sigmund posited that the local erosion rate of the surface is 
proportional to the local atom emission rate resulting from the atomic collision cascade, 
and that the emission rate at a point on the surface is proportional to the nuclear energy 
deposition density at that point resulting from collision cascades from the ions impinging at 
all points. Sigmund subsequently [35j recognized the destabilizing influence of the curvature- 
dependence of the sputter yield (atoms out per incident ion) by modeling the nuclear energy 
deposition density as taking the form of Gaussian ellipsoids beneath the surface and showing 
that, as a consequence, concave regions of the surface receive more energy and thereby erode 
more rapidly than do convex regions 57l |. 

The origin of the characteristic length scale of the self-organized patterns was identified 
by Bradley and Harper (BH) 4|], who recognized that Sigmund's destabilization mechanism 
is opposed by surface diffusion, which operates so as to return the surface to flatness 58|]. Ex- 
panding Sigmund's Gaussian ellipsoid response in powers of derivatives of the surface height 
h{x, y, t) and superposing classical Mullins-Herring [36|, l37| surface diffusion, BH derived 
a linear partial differential equation (PDE) {4] that describes the evolution of the surface 
height on scales much larger than the characteristic length scales of Sigmund's Gaussian 
response: 

-I + {SAx + Sydyy " BV^}h , (l) 



dh 
'di 



where I{h) is the vertical erosion rate of a flat surface, Sx{h) and Sy{h) are the curvature 
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FIG. 1: (a) Plot of sputter yield curve I{0), normalized by /(O) (b,c) Plots of Sx{0) and Sy{9), 
normalized by |5x(0)| = 152^(0)1. The parameters used are: a = 1.5 nm, a = 0.9 nm, /u = 0.5 nm. 



coefficients, b is the surface slope, and i? is a material parameter describing relaxation and 
containing the surface diffusivity and the surface free energy. The coefficients I,Sx, and 
Sy are expressed in terms of Sigmund's Gaussian and depend on 6' = tan~^(6), the angle 
between the beam direction, henceforth denoted as —z, and the local normal to the surface 
n (0 < 6 < 7r/2). For nonzero 6, we denote by x the axis perpendicular to z in the n — z 
plane. Bradley-Harper's linear stability analysis yields unstable modes whenever Sx or Sy is 
negative, whose characteristic length scale arises from a balance between the destabilizing 
effect of the second derivatives dxx-, dyy and the stabilizing effect of the surface diffusion term 
V^. The behavior of Sx{0) and Sy{9) for characteristic parameter values are shown in Fig. 
1. The Bradley-Harper analysis gives rise to the following predictions: (i) Below a crossover 
angle 6 cross, Sx < Sy < 0, implying a faster growth rate for parallel mode (wave vector 



parallel to projected ion beam direction a 
(wave vector along y) surface modulations 



ong the surface) than for perpendicular mode 



591 ]: (ii) Sy < for all 6, implying instability to 



perpendicular modes at all incidence angles. For 9 > Ocross the perpendicular modes are 



the fastest to grow with dominant wavelength ^JSn'^B /{—Sy). The generalization of the BH 
analysis to the nonlinear regime, which is required to account for the observed saturation of 
ripple amphtude and the emergence of more complicated patterns (e.g. hexagons, dots, pits) 
was carried out by Cuerno and coworkers [s], [2l| who expanded Sigmund's Gaussian ellipsoid 
model to higher order in surface height derivatives, resulting in a Kuramoto-Sivashinsky type 
equation [38| for the surface evolution. 

There is growing evidence that although the Bradley-Harper predictions explain some 
features of experiments (e.g. the temperature dependence of the wavelength of the ripples 
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[l3|), there are also some glaring inconsistencies. This is clearly demonstrated in, e.g., the 
recent work of Ziberi et al. [31], who found a "window of stability" for Si surfaces at room 
temperature bombarded by ~ 1 — 2 keV noble gas ions at an intermediate range of angles 
di < 9 < $2, where 6i ^ 30°, and 62 ~ 60°. Moreover, Ziberi et al. demonstrated that 
when bombarded by some noble ions (Ne"*"), a fiat surface remains stable at all angles. In 
addition to the experimental inconsistencies with BH prediction (ii), there have also been 



recent experiments 
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391] and atomistic simulations 40|, |4l| that have measured the shape 



change of a smooth solid surface in the vicinity of an impingement by a single energetic 
monatomic ion or cluster ion. These studies show significant deviations from the predictions 
of Sigmund's ellipsoidal Gaussian form. For example, the molecular dynamics studies of 
Feix et al. [27j] indicate that for 5 keV Cu^ bombardment of Cu crystals, the collision 
cascade intensity along the surface has a maximum along an annulus some distance from 
the impact point and its spatial decay is better characterized by an exponential rather 
than by a Gaussian function. In this case Feix et al. still found linear instability of a flat 



surface. Moreover, in many cases 
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401 ]. including low energy (0.5 keV) bombardment 



of an amorphous silicon surface 4l|], the response of the surface is the formation of craters 
with rims. This type of response, involving the accumulation of matter at some locations, is 
in clear contradiction to the purely erosive response predicted by Sigmund's model using a 



Gaussian ellipsoid co 



lision cascade. The occurrence of craters with rims has been attributed 



to thermal spikes 40|] or to ion-stimulated surface mass transport 41]. 

These observations raise the interesting question of how robust are the predictions of 
BH to the precise shape of the local response to an ion impact. Indeed, the most general 
evolution equation based on the accumulation of local responses to ion impacts is 33] 

(9/i(x,t) 



dt 



J dx' Jio„(x') A[x - x, h^{x, t), hy{x, t), h^x{^, t), hyy{x, t), h^yix, t), 



(2) 



where x = {x,y), Jion{^') is the ion fiux at x', subscripts x and y denote partial derivatives, 
and the kernel A[x — x', . • • ], representing the change in height at x due to an ion impact at 
x', is expected to decay smoothly to zero at large distances |x — x'|. This equation is more 
general than that assumed by Sigmund because the kernel A can have any shape whatsoever, 
and can depend on the complete local geometry of the surface. 

In this paper we explore whether a more general physically motivated surface response can 
change the predictions for linear stability from those of Bradley and Harper. Our purpose 
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here is not to perform quantitative comparison between theory and specific experiments, 
but rather to determine how robust the predictions of the Bradley-Harper theory are with 
respect to modifications of the ion impact function A. We demonstrate that, whereas the 
fundamental prediction concerning the instability of flat surface to uniform ion irradiation 
results from a wide class of response functions including Gaussian and non-Gaussian dis- 
tributions - thus explaining the applicability of Bradley-Harper theory for wide range of 
systems - there are certain classes of modification that have a dramatic effect. Notably, 
these modifications render the fiat surface stable - in contradiction to the classical theory - 
while imperceptibly affecting the yield curve I{b). 

The paper is organized as follows: In section II we extend the BH approach - of deriv- 
ing from the microscopic response function the coefficients S^i^b), Sy{b) in Eq. ([T]) - to a 
broad class of purely erosive surface response functions, of which the Gaussian ellipsoid is a 
particular example and the response ofFeix.ol, 0, sanotK. example, we sKow that 
the BH prediction of linear surface instability for all incident beam angles is unchanged. 
Hence any purely erosive surface response within this broad class is contradicted by exper- 
iments. In the remainder of the paper we explore possible physical mechanisms that could 
resolve this condundrum. In section HI, we demonstrate that a surface response that is not 
purely erosive, but rather consists of the formation of a crater surrounded by a rim, does 
allow linear stability for some range of incidence angles. In section IV, we demonstrate that 
impact-induced "downhill" surface currents, such as those recently found in MD simulation 
of C and Si surfaces bombarded by low energy (~ 250 keV) ions 42], can also yield linear 
stability for some range of beam angles. There are thus multiple physical mechanisms that 
could explain the experiments, and the essential question is to determine which effect is 
dominant. Identifying the dominant physical mechanism for linear (in) stability is critical 
to having a reliable nonlinear theory for pattern formation. In section V, we discuss how 
experiments might distinguish the competing theories. In particular we argue for a careful 
analysis of experiments near the observed critical angle at which a fiat surface becomes 
stable. 
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II. BRADLEY-HARPER THEORY REVISITED 



The Sigmund theory of sputtering [SJ] posits that the local erosion of the surface in 
Eq. ([2]), A[x — x']/vTT~P^ with b the local surface slope, is proportional to the local atom 
emission rate resulting from the nuclear collision cascade, which itself is proportional to 
the nuclear energy deposition density at (j£,h{x.)) from an ion impinging at (x', /i(x')). To 



demonstrate the source of an instability 351], Sigmund modeled the collision cascade as a 
Gaussian ellipsoid. Bradley and Harper's subsequent expansion of Sigmund's Gaussian el- 
lipsoid collision cascade model, combined with smoothening by fourth-order Mullins- Herring 
surface diffusion, leads to Eq. ([T]). 

To examine the consequences of forms of the erosive response that are more general than 
Gaussian ellipsoids, we assume: 

A[x-x', .••] = Ah{r,z) 

= -Ae-^^'^)-^^^) , (3) 



where r = i^xM-^, z = h{x,y), and A is a length that depends on parameters such as 
ion energy and ion and target mass. The first equality in Eq. assumes radial symmetry 
about the ion track and no explicit dependence on the surface slope and curvature, with the 
kernel depending only on r and z. The second equality assumes separation of the variables 
r and z. In Eq. ([3]) the ion is assumed to penetrate the surface at (r, z) = (0, 0). 
Sigmund's Gaussian ellipsoid response is a particular case of Eq. ([3]), with 

f{z) = ^,{z-ay ; 9ir) = ^,r\ (4) 

where a is the average penetration depth of the ion, and a, ^ are lengths characterizing the 
ranges of response in directions parallel and perpendicular to z, respectively. 

Following Bradley-Harper, we substitute in Eq. ([2]) the response form ([3]) and add a 
relaxation mechanism to the surface dynamics associated with Herring- Mullins surface dif- 
fusion: 

^ = -BV^h - a Hdy Hdx e-9(V^)-m^^y)) ^ (5) 

where a = AJion and the materials parameter B is given by i? = '-/Q'^DC /{kBT). Here C, 
D, and fl are the concentration, diffusivity, and volume, respectively of the surface- diffusing 
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species; 7 is the surface free energy, fc^ is Boltzmann's constant and T is the absolute 
temperature. 

To study evolution of surface morphology in the limit that the surface height h{x, y, t) 
varies on scales much larger than the ion penetration depth, we consider perturbations about 
a planar surface (x, y,h = bx), so that 

1 1 

h{x, y) = bx + -KxX^ + -hyyy"^ + h^yxy H , 

and expand e'^^^^^'"^^^ to obtain 

exp[-f{hix,y))] ^ e-^(^^)[l - f'ibx)[h,,x^ + hyyy^ + KyXy]] . (6) 

With the expansion ([6]), the integral equation ([5]) is readily transformed into the PDE ([1]) 
with the coefficients: 
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I{b) = a dy dx e-^''(^'^) 

J —00 J —00 
POO POO 

Sy{b) = a dy dx e-f"'^''^y^f'{bx)y 



'OO J —00 
00 poo 



S^{b) = a dy dx e-P^^''^y^f{bx)x^ (7) 



where pb{x,y) = g{\/ x^ + y^) + f{bx). 

The question now is how various choices of /(r) and g{z) can change /(&), Sx{b) and Sy{b). 
We are primarily interested in the slope dependence in Sy{b), because in the Bradley-Harper 
theory Sy{b) < for all slopes b. Our question is whether any choice of f{z),g{r) can stabilize 
the surface against perpendicular modes {Sy > 0) for some range of b while not significantly 
affecting the shape of the yield curve. The latter requirement is especially significant because 
the yield curve predicted by the Sigmund response function agrees quahtatively with that 
measured on many materials - at least for non-grazing incidence {3]. 

All of our analysis proceeds with the same methodology: the integral for Sy{b) in Eq. 
([7]) is dominated by contributions near the minimum of pt which we call {xminiUmin}- This 
is because the size of the region where energy is deposited (of order the penetration depth 
a) is much smaller than the characteristic length scale over which the surface shape varies. 
The minima of pb satisfy the equations 



Urain '(12 i 2 



9'{\ hlrin + yLn) = 0; (8) 



2 ' 
min 



''■'^Ln + yLn)+M'{hXrmn) = 0. (9) 



Depending on the functional forms of g and / there are two possible types of solutions to 
these equations. 

(a) Vmin = 0, ±g'{Xmin) + bf'{bXmin) = (10) 

(b) 9'i\/xLn + yLn) = 0, &/'(&a;™n) = , (11) 

where the ± signs in (ITUi) correspond to Xmin > 0, Xmm < 0, respectively. Once the locations 
of the minima are determined, we can expand 

Xmin, Vmin ) H (fc + b f ) 

, ( _ V _ _l_ ~ Vmin) 

~r \X Xjninjyy ymin ) Qxy ~r ^ Qyy 

= P* + A{X - Xminf + C{X - Xmin){y " 1/mm) + B{y - Uminf (12) 

where the second equality defines A, i?, C, p* . This expansion can then be used to evaluate 
the integral. 

We now proceed to use this methodology to establish the conclusion that 5*^ < is ex- 
tremely robust. For any kernel of the form considered here a perpendicular mode instability 
always exists for all slopes b. The characteristic behavior of the coefficient is more fickle. 
Obviously, for 6^0, S^ib) / Sy{b) 1, and therefore S^ib) is necessarily negative for small 
enough slopes b. However, Bradley-Harper's observation, that Gaussian ellipsoids imply 
Sx < Sy < Q for 6^1, does depend on the exact shape of the response function. This can 
be readily verified by considering Sigmund's response Eq. (jlj) with a < a. Hence, we will 
focus our analysis on the robust properties of the linear dynamics, associated with the sign 
of Sy, and will not further discuss in this section. 



A. The shape of the energy distribution does not quaUtatively affect stabihty 

We begin by considering changes in only the shape of the energy distribution: namely 
we consider f{z),g{r) that keep the position of maximum energy deposition at a single 
point (the average stopping point of the ion), though we vary the shape of the distribution. 
We thus assume that the function f{z) has a minimum at z = a whereas g{r) increases 
monotonically from r = 0. 

Under these assumptions, the minimum of pb{x,y) must be of type (fTOj) . Moreover, 
because the minimum of g{r) along the x axis occurs at x = and the minimum of f{bx) 
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occurs at X = a/h > 0, then Xmin, determined by g\xmin) + bf'ihxmin) = must be in the 
interval < Xmin < a/h, such that f'ihxmin) < 0. The expansion of ph in equation f|T2l) leads 
to the coefficients A = {g" + 6^/")/2, B = g'/2\xmin\ and (7 = 0, where all derivatives are 
taken at Xmm- Hence the integral is approximately 

/oo POO 
dy / dx e-^*-^^^~^--^'-^y'nhXrmn)y'. (13) 
-co J —oo 

Because f'{bxmin) < the integral ( |T3l) is necessarily negative for all b. This demonstrates 
that the experimentally observed stability of a sputtered surface to perpendicular mode 
ripples is not a consequence of the shape of the energy distribution. 



B. Toroidal energy distributions do not qualitatively affect stability 

Another possible modification of the energy distribution is for the maximum energy 
deposition to occur away from the ion trajectory. Indeed, Feix et al.'s recent simulations 



of Cu crystals bombarded by 5 keV Cu+ ions |27| have demonstrated energy distributions 
with a maximum along an annulus surrounding the ion trajectory. Such a response is thus 
characterized by a g{r) with a minimum at Vmin = ''"o > 0. 

Consider the sign of Sy{b) under these circumstances. There are now two different regimes, 
depending on the slope. When the slope is small, such that a/b > tq, the minimum must 
be of type (a), Eq. ffTOj) . Type (b) (Eq. f|TT]) ) is excluded because if fibx^in) = then we 
must have Xmin = a/b. But then the equation g'{r) = cannot be satisfied: this equation 
imphes that + = Tq, which cannot be obeyed for any ymin- In contrast, when the 
slope is large, so that a/b < tq, the minima are of type (b). 

Let us ffist consider the regime of small slope. Here the analysis proceeds as above with 
the same A, B, C defined in f|T2l) . As before the sign of the integral hinges on the value of 
f'ibxmin) = —g'{.Xmin)/b. Bccausc wc are assuming that the minimum of f{bx) occurs at 
X = a/b which is larger than the minimum assumed by g{r) along the x-axis, at x = tq, Eq. 
(|TOl) implies that f'{bxmin) < 0. Hence we arrive at the conclusion that in the small slope 
regime Sy < 0: the linear instability survives. 

The second regime, where b/a < vq, is more subtle, with two minima being of type (b) 
(Eq. [TT]) . Assuming the minimum of g{r) occurs at Tq, and the minimum of f{z) occurs at 
a, in this case we have that {xmin, Umin) — ('^/^? ±A/rQ — (a/by). The value of Sy{b) is given 
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by the sum of the contributions to the integral centered around each of these two minima. 
For these minima the values of A, B, C are given by 

where g^' is evaluated at r = ro and /" is evaluated oX z = a. We now must evaluate 

/oo nOQ 
dy I dx e-''*-^^^-^— )'--^(^"^-")'-^*(^-^™-)(^-^-")f (15) 
_^ -oo J —oo 

The exponential in the integrals are best dealt with by completing the square, so that they 
become 

/oo poo —p*~A\ (x — Xmin) + ^-^(V — V^- ) I 

dy dxe \ ) e-^y-yt.SiB-^c^fl^^)fi^^^)y\ 

_i_ -oo J —oo 

(16) 

Now the second exponential decays with y varying away from y^^^ because B > ^ for 
any b ^ 0. If we now change variables to x = x - x™„ + - Vmin) and y = y - y^i^ we 
obtain 

/oo fOO 
dy cixe--^^^"^^(^-(^")V4^)/'[6(x_ + x-^^)](l/±„^ (17) 
_^ -oo J —oo 

Because now f'ihxmin) = 0, evaluation of the integrals to leading order requires expansion 
of the terms f'[b{xmin + x — ^y)] around a = hx^in- With this we get the following 
approximation to the integral: 

/oo /"OO / \ 

dy / di e-^^'-fiB-i^C^?I.A) _ ^ . . . ^ 

_i_ oo J — oo V / 

(18) 

The contribution of the two integrals is identical, and sums up to: 

S^{h)^-na)g"{T,)'^V, (19) 

rlA 

where we have substituted the formula for (fMl). have used hxmir, = a, and where 



' 5 

fOO POO 



/OO POO 
CO J — oo 

The RHS of Eq. (1191) is negative definite, i.e. Sy{b) is negative for all values of b. Hence 
response functions of the form of Eq. generally cause a perpendicular mode instability 
for any incidence angle. The qualitative conclusions of the original Bradley-Harper analysis 
concerning the instability of perpendicular surface modulations at any beam angle are thus 
very robust. 
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III. EFFECTS OF MASS REDISTRIBUTION 



The analysis of the previous section demonstrates that a broad class of purely erosive 
response functions gives rise to linear instability for all beam angles. However, there have 
been several recent studies suggesting that the surface response is not purely erosive. These 
studies demonstrate that after ion impact, a crater forms around the impact point of the 
penetrating ion, surrounded by rims elevated from the original surface 18|, |39|, |40|, |41| . 
This behavior, where Ah > in the rim, is completely different from the erosive response 
functions described above. We investigate whether such response functions can cause the 
stability of a flat surface. 

To carry out this analysis, we introduce a natural generalization of the family of response 
functions ([3]): 

A/i(r, z) = -Y^ A^.g-f^M-/^ W (20) 

where gj{r),fj{z) are localized functions as discussed in the previous section, but the co- 
efficients Aj can be negative or positive. In particular, negative Aj corresponds to mass 
deposition associated with ion impact and can give rise to formation of rims. A particularly 
simple form of a response function is the sum of two Gaussian ellipsoids: 

This response function has eight free parameters (including A and (3), all of which are con- 
strained to be positive. Unlike the original Sigmund model, the free parameters here are not 
directly connected to a microscopic picture. Because our intent is to understand whether 
small deviations from Sigmund's response function can change the stability characteristics 
of the surface, we will consider the case with /5 -C 1, and think of ai, yUi, cti as corresponding 
essentially to the original Sigmund parameters. The parameters a2,/i2,cr2 describe charac- 
teristics of the mass redistribution. 

With the model so defined, we can evaluate the yield curve I{b) as well as S^{b), Sy{b), 
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obtaining 



i=l,2 



j=l,2 V t^i+^i) 



Sy{h) = -2^J_ Y ^^^^^^^n/ ' (24) 



j=l,2 



where we used the notation Ai = A,A2 = —13 A. 

We now want to use this result to address the following question: is there a regime of 
parameter space where the stability characteristics of the surface are qualitatively different 
from the predictions of Bradley and Harper, but for which the yield curve is experimentally 
indistinguishable from that predicted by the Sigmund response? Indeed, we have found 
multiple regions of parameter space where this occurs. This can be demonstrated simply 
and analytically by expanding equations ( 122]|23|I24I) in the regime of small slopes, where 
5*^. Sy. We find that, as 6 0, 

I{b) ^ 27rJi„„A[/i2e-'^?/(2'^?) - /?/i2e-'^2/(2-i)] (25) 
Sy{b) ^ Sxib) ^ -27rJ_A[^e-"?/(^-?) - p^e-y^'^'^^] , (26) 

Here we see that for small slopes, and Sy can have either sign, depending on the relative 
magnitudes of the terms ^^i^e~"i'''-^'^i'' and /3^2^e~'^2/(2o-2)^ jf ^j^g second term dominates the 
first then and 5*^ are positive at small b and the surface is stable to all perturbations. Can 
stability be achieved without significantly affecting /(&)? Obviously, this will be the case 
if /i^e~"i/*^^'^i) ^ (3 file'"'^'^'^'^^) , Letting Zi = /ife"'^?/*^^'^?\ satisfaction of the two conditions 
amounts to finding parameters where (i) Zaaifil/af < l3Z2a2n'i/ ai while (ii) Z1/Z2 ^ (3. 
We also would like (3 to remain small. Such a parameter regime clearly exists and merely 
constrains the scale and the geometry of the mass redistribution region. 

To demonstrate this explicitly. Fig. [2] shows the behavior of I{b), S^^b), Sy{b), where we 
have used the same parameters for ai = 1.5, ai = 0.9, /ii = 0.5 as used for the 'normal' 
Bradley-Harper stability characteristics shown in Fig. [H with the additional parameters 
f3 = 0.03, 02 = 0.5nm, (T2 = 0.5 nm, and /i2 = 1 nm. For these parameters Zi = 0.06nm^ 
and Z2 = 0.6nm^, thus Z1/Z2 ^ (3, whereas aifif/af = 0.46 nm and 02/12/0"! = 2 nm. 
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FIG. 2: Normalized yield curve and BH coefficients Sx,Sy for two sets of parameters of the two- 
Gaussians model, Eq. (j2ip . The "Sigmund parameters" a,a,^ are taken as in Fig. [H and the 
same normalization factors are used. The new parameters are: (a) a = 0.03; a2 = 0.5nm,o"2 = 0.5 
nm, and /i2 = 1 nm and (b) a = 0.03, 02 = 0.9nm, cr2 = 0.2 nm and /i2 = 1-5 nm. 



We therefore satisfy both constraints (i) and (ii) listed above. Indeed, the top row of Fig. 
2 shows a stable region of parameter space for small slopes in both Sx and Sy, while the 
qualitative shape of the yield curve is unchanged. We have also found regions of parameter 
space where the two conditions derived above are not met, hence a flat surface is unstable 
at small b, but there is still a window of stability at higher slopes, as shown in the bottom 
row of Fig. 2. 

The results of this section demonstrate a very significant conclusion: that small changes 
in the shape of the surface response of a single ion can completely change the stability 
characteristics of a flat surface from those predicted by Bradley and Harper, but yet not 
lead to any significant modification to the measured yield curve. Further analysis along 
this line requires a microscopic theory for the non-erosive processes, or detailed atomistic 
simulations from which effective parameters such as /3, a', a', and fi' can be determined. 
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IV. INDUCED SURFACE CURRENTS 



In the previous sections we considered a surface response that does not depend exphcitly 
on the incidence angle and is fully characterized by considering normal incidence (6 = 0). 
Namely, the response at a point {x,y, h{x,y)) depends only on the projections of the vec- 
tor that connects {x,y, h{x,y)) to the average ion stopping point (0,0, —a), in directions 
parallel and perpendicular to the beam direction z. Thus, the dependence of the coeffi- 
cients I{b),Sx{b) and Sy{b) in Eq. ([1]) on the angle 9 = tan~^(6) is imphcit and purely 
geometrical, stemming from the fact that the distribution of values of these projections 
+ a| , ^/x'^ + 1/2, respectively) over all surface points depends on the slope b. 

It is possible, however, that the response of a surface point to ion impact depends exphc- 
itly on the incidence angle. Such behavior was reported by Moseler et al. [42], who used 
molecular dynamics to study the ion-enhanced smoothening of diamond-like carbon sur- 
faces bombarded by low energy (30-150 eV) carbon ions. These authors simulated surfaces 
tilted at angles up to 20° and observed transient surface currents with components along 
the projection of the ion beam direction onto the surface, resulting in net displacements 
along the surface of magnitude proportional to the incidence angle. Their analysis of this 
effect, neglecting densification and sputter erosion, and focusing on beam angles near normal 
incidence, resulted in an isotropic diffusion-like equation for the surface height: 

dh/dt = uV^h , (27) 

where u is positive and consequently stabilizing (c.f. Eq. (1)). Moseler et al. did not pursue 
the beam angle dependence of u. As is the case for the erosion coefficients in Eq. ([T]), we 
expect this smoothening effect to become anisotropic away from normal incidence, yielding 
two different coefficients z/^(6), i^y{b). 

Previously, Carter and Vishnyakov proposed a similar smoothening term to explain 
the absence of linear instability on silicon bombarded with 10-40 keV Xe"*" at incidence 
angles between and 45°. They proposed a mechanism whereby forward recoils move, 
on average, parallel to the ion beam before coming to rest. They retained the projection 
along the surface, which may be interpreted as a consequence of the incompressibility of 
the solid: the surplus density injected into the solid subsequently "pops up" to the surface 
along, on average, the shortest path. Specifically, for an ion flux of magnitude Jjo„ in a 
plane perpendicular to the ion beam, the number of ion impingements per unit area of 
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surface is JionCos{9), where 9 is the local angle of incidence. The induced current per ion 
projected along the surface varies as sin(^), resulting in a surface current Jx proportional 
to Jjon sin(6') cos(^), or Jjo„sin(2^). This surface current has the same stabilizing effect on 
parallel mode instabilities as that identified in the simulations of Moseler et ai, Eq. flTTI) . 
but with Ux oc cos(2^^). Carter and Vishnyakov did not consider Uy. 

In principle, the low-energy mechanism of Moseler et al. differs from the high-energy 
Carter- Vishnyakov mechanism: in the former case, the projected range is ~ 1 nm and 
true surface transport is observed; in the latter case, the projected range is greater than 
10 nm, volume transport is induced, and it is the component parallel to the surface that 
results in the smoothening effect. However, in both cases an explicit dependence on angle 
of incidence is apparent, and phenomenologically they appear virtually indistinguishable. 
In both mechanisms the average net effect of each ion impact is a displacement along the 
surface that is proportional to 9 for small 9 and should saturate at large 9, as does sin(^). 
In all cases the ion impingement rate per unit area of actual surface goes as cos(^). Their 
combination should result in an induced "downhill" surface current that approaches zero 
near normal and grazing incidence and displays a maximum in the vicinity of 45°. 

To understand the implications of Eq. fl27|) for linear stability, it is essential to establish 
the dependence on incidence angle of both coefficients i'x{9), i^y{9) for parallel and perpen- 
dicular modes, respectively. To this end we consider a simple model in the spirit of those 
discussed above. The geometry of the previous sections is assumed, where an ion flux Jjon 
impinges in the —z direction on a surface slightly perturbed from the plane h{x, y) = bx, 
and 9 is the angle between the local normal to the surface and the z axis. We assume that 
the component of ion momentum parallel to the surface causes the displacement of surface 
target atoms a distance along the surface proportional to sin(^). The contribution of the 
induced surface current Jg = {Jx,Jy) to dh{x,y,t)/dt is —V • Jg, where V = {dx,dy). In 
order to evaluate Jg let us assume first that the surface is exactly described by h{x, y) = bx, 
where b = tan{9). In this case Jy = 0, and with a momentum component parallel to the 
surface proportional to sin{9), we obtain Jx oc — J^on cos(6') sin(6'), where JionCos{9) is the 

rate of ion impingement per unit surface area. This behavior is consistent with the results 

' ii 

of the MD simulations of Moseler et al. (j42?]). In order to write the induced surface flux 
for a general surface, represented by the equation z = h{x,y), we must express Jx and Jy 
in terms of Vh. The angle 9 satisfies the relation cos(6') = n ■ z = l/v'lV/iP + 1, where 
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FIG. 3: Normalized coefficients Sx^i^x^Sy, and Uy, comparing the effect of surface-induced currents 
{u, Eq. [27|) to erosion from Gaussian ellipsoids (S*, Eq. [1]). The relative magnitude of v/S at 
normal incidence varies with conditions and is chosen arbitrarily here for illustrative purposes. 



h = [—dh/dx, —dh/dy, 1]/ a/| V/ip + 1 is the unit vector normal to the surface. Let us de- 
note by the angle between x axis and the direction within the x — y plane of maximal 
increase in surface elevation at {x, y): (j) = tan~^ M/lf • "^^^ fluxes J^., Jy are then given by 
Jx oc — sin(26') cos(0) and Jy oc — sin(26') sin(0). 

Because our analysis in this paper is restricted to linear dynamics of the surface, we 
expand V ■ J to linear order in deviations of h from the flat surface h = bx {b ^ 0). 
Algebraic manipulation yields the relations: 

cos(0) ~ 1 , sin(0) ^ b — , 

oy 

b f)h b 1 — /)2 

cos(^)- (1 + 6^/^(1-^-) , sin(2^)^^[l + ^^^aV5a:], (28) 

and the linear contributions i^xib), z/y(6)from the surface induced currents to the coefficients 
Sx{b), Sy{b), respectively, in Eq. ([T]) is: 



Ux{b) oc 
Uy{b) oc 



1-b^ 

1+62)2 
1 



(29) 
(30) 



1 + 62 

The expression for Ux in Eq. (I30l) is equivalent to the expression derived by Carter and 
Vishnyakov [60|] . Notably, the mechanism described by Eq. (1271) corresponds to a conserved 
surface current and thus does not have any effect on the yield curve I{b). The effect of 
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induced surface currents on the stability is evident in Fig. [31 The effect stabihzes both modes 
from normal incidence up to incidence angles of 45°, whereupon it becomes a destabilizing 
influence on only the longitudinal mode. The magnitudes of and Uy must equal each 
other at normal incidence, but their relationship to the magnitudes of Sx and Sy depends 
on the relative strengths of the mechanisms. If the induced surface current mechanism is 
sufficiently strong, as illustrated in Fig. [3l then starting with normal incidence and going to 
increasing angles, one should observe a regime of absolute stability; the dominance of parallel 
modes; and the dominance of perpendicular modes. For further insight, it is essential to 
estimate the strength of the induced surface current and how it depends on materials and 
ion beam parameters, e.g. by methods such as atomistic simulations. 



V. EXPERIMENTAL SIGNATURES 



We have described several mechanisms by which surface dynamics of the form ([T]) can 
account for regions of ion beam angle where a flat surface can be stable or unstable. The 
mechanisms suggested in the previous two sections provide some scenarios leading to mod- 
iflcations of the Bradley-Harper coefficients in Eq. ([1]) and thereby causing stability of 
the bombarded surface at various ranges of angles; there are also potentially other such 
mechanisms. 

The critical question now is to determine which of the potential physical effects is operat- 
ing in experiments; the answer to this question almost certainly depends on the material, the 
ion mass and energy, etc. Beyond the linear stability analysis itself, this issue is of central 
importance for developing a quantitative nonlinear theory of pattern formation; it is well 



known [38| that accurately identifying the linear dispersion relation is critical for deriving a 
nonlinear theory which can predict the fully developed pattern. 

How can experiments discern the dominant linear (in) stability mechanism? Here we 
present one method for ruling out some of the possibilities: in particular we point out 
the relevance of the stability-instability transition not only as an interesting dynamical 
phenomenon, but as a conceptual tool to gain valuable information on the general character 
of the dynamics of ion sputtered surfaces further away from the transition. 

In general, the linear stability analyses discussed in this paper result in a dispersion 
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FIG. 4: Schematic plots depicting the transition between stable and unstable surface dynamics 
for three dispersion relations, (a) Left column: generalized Bradley-Harper, Eq. (I3ip . where 
the transition occurs at S'^^'* = with diverging wavelength, (b) Middle column: with Facsko 
nonlocal "damping term", transition occurs at S^^'* < with finite wavelength, (c) Right column: 
with Asaro-Tiller nonlocal elastic energy mechanism, transition occurs at S"^^'* > with finite 
wavelength. 



relation of the form: 

R, ^ Re (cuq) = -S^ffql - Sfql - B^^ql - Byyq^ - B^yq^q^ + ■ • • , (31) 



6l| 2l| which describes the growth rate of a Fourier mode: 

V,..W=^.,.«(0)e^^^^^+''"^^+"''*- (32) 

In equation ([31]) we have lumped the two quadratic contributions into S^-^J^ = Sx,y — i^x,y We 
focus on the transition between stable and unstable perpendicular (parallel) modes described 
by Eq. fl?I]) as Sy-^^^ (S^-^^) changes sign. This is depicted in the left column of Fig. (4). 
Here we assume for simplicity that the only parameters in Eq. (1311) that change appreciably 
with the beam angle are S^'^-^ , Sy^^ . 
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The most important feature of this schematic plot is that it predicts divergence of the 
pattern wavelength upon reaching the transition to stable surface dynamics. To see this more 
clearly, notice that a condition for the existence of linearly unstable modes is that ma.x{Rq), 
the maximal value of Rg over all wave vectors q = [q^, Qy) is positive. Assuming a smooth 
dependence of all coefficients on the beam angle, a transition between stable and unstable 
surface dynamics corresponds to a beam angle for which max(i?g) = 0. For simplicity, let 



us assume that max(i?g) is achieved for q = {qmax,0). Then: qmax = y —Sx /^B^x and 
max{Rq) = R(q^,,fl) = -Sff/AB^^ = 0, implying S^^f{9) ^ S* = and hence Qmax at 



the transition. 

A diverging length scale is a strong characteristic signature of the stability-instability 
transition, and it is thus natural to ask whether this prediction is valid if other physical 
processes, not accounted for in this paper, influence the surface dynamics and thus modify 
the dispersion relation (I3T!) . We argue that this divergence is expected as long as the 
following assumptions are satisfied: 

1. The beam-angle dependence of all coefficients in the equation is smooth. 

2. The linear dynamics is analytic, ruling out terms like \Vh\ 

3. The dynamics is first order in time. 

4. Linear surface dynamics is local - namely, it can be described by a partial differential 
equation (PDE). 

Assumption 1 is required because, as can be seen easily from Fig. (4a), a discontinuous 
"jump" between negative and positive values of S^^^^Sy^^ at some beam angle 9* may 
yield a transition to stable dynamics at |g| > 0. Physically, a discontinuous change of 
parameters is associated with abrupt changes in material properties, such as amorphization 
of a crystalline surface. For such a scenario to be associated with a smooth change of the 
beam angle is sufficiently unlikely as to be a rare occurrence. Assumption 2 is required in 
order to make a linear stability analysis meaningful. If this assumption is violated then the 
early stage surface dynamics of an initially flat surface is not described by the dynamics of 
independently evolving Fourier modes (1321) . Assumption 3 is expected to hold as long as 
inertia is neglected. Assumptions 3 and 4 together imply that the the ampliflcation rate 
Rq, which is the real part of the complex eigen- frequency a;^, contains only even positive 
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powers of q. Namely, local processes, by which a change of surface height is related to the 
variation of erosion or flux rates between a surface point and its nearest neighbors can be 
described by spatial derivatives of the function h{x,y,t). In a dynamics that is first order 
in time, the eigen-frequency Ug in Eq. fl52]) thus equals a polynomial in q, where all spatial 
derivatives with odd order (i.e. Vh, V^h) have imaginary coefficients, and thus do not 
contribute to the amplification rate Rq = Re (ujq). Notice also that the locality assumption 
rules out the existence of a constant term (i.e. oc g°) in fl32l) . This is a consequence of the 
invariance h ^ h + const.. Therefore, a term oc h{x,y,t) (i.e. without spatial derivatives) 
can appear in the surface dynamics only as a combination respecting this invariance such as 
h{x,y,t) — h{t), where h{t) = f (ix/i(x, t), and thus must be associated with some nonlocal 
processes. 

Thus, under these general assumptions (and neglecting the possibility that spatial deriva- 
tives of order 6 or higher are dominant in the dynamics), the amplification rate Rg satisfies 
Eq. ( !3T|) . the stability- instability transition is depicted by Fig. (4a), and the characteristic 
wavelength at the transition is predicted to diverge. 



Recently, Ziberi [4^ and George [45|] have measured the pattern wavelength at several 
values of beam angles near the transition to the stable region in silicon irradiated by no- 
ble gas ions at temperatures where the surface should be amorphous and isotropic. The 
measurements indicate that the wavelength at the transition remains finite, and may thus 
be a strong indication that one of the above assumptions is violated. Anticipating that 
assumptions 1-3 are still valid, we will discuss here two nonlocal terms, whose introduction 
may render the wavelength at the transition finite. 

A. Facsko "damping" term 

First, let us consider the effect of including a linear term, K[h{x, y, t) — h{t)] with h{t) = 
J (ix/i(x, t), in the surface dynamics, Eq. ([1]). Such a term was recently introduced by 



Facsko et. al. [23| as a possible way to obtain long range ordered patterns observed in the 
fully nonlinear regime. The term is suggested to be a placeholder for a model of redeposition. 
With such a term, a constant K is added to the right hand side of the dispersion relation 
pril . This is consistent with the dispersion relation measured by Brown and Erlebacher 
291] on Si(lll) at temperatures where it should remain crystalline, with singular surface 
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energetics (making the validity of assumption (2) questionable). The effect of this term 
on the transition between stable and unstable dynamics is depicted in the middle column 
of Fig. (4), where it is demonstrated that the characteristic wavelength does not diverge 
at the transition, as can be obtained from the following analysis: Again, for simplicity we 
assume that max{Rq) is achieved for q = {qmax,0)- Here again q^ax = \J —St/^ /2Bxx but 
max(/?g) = R(g^,,,o) = K -Sl^^ /AB^^ = 0, implying S^^f{9) = S* < and hence \qmax\ > 
at the transition. 



B. Asaro-Tiller mechanism 



The Asaro-Tiller elastic energy driven mechanism 46|, |47| gives rise to instability of solid 



surfaces under biaxial in-plane stress. Biaxia 



in the bombarded solid 
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compressive stresses are known to develop 



54| . and this effect could be important in 



the surface dynamics. Assuming a sinusoidal modulation of a free surface under biaxial 
compressive stress, the tangential stress increases at the troughs (compression) and decreases 
at the peaks (dilation) by an amount proportional to the wavenumber of the modulation 
and to the applied stress uo in the solid. This increases the chemical potential at the 
troughs compared to the peaks and drives a surface current from the troughs to the peaks 
that further amplifies the modulation, thus leading to instability. Including this effect in 
the surface dynamics gives rise to a term oc M|gp on the RHS of Eq. ( pTll [sS^, where 
M oc (jg. This term does not stem from local effects but rather from nonlocal effects 
associated with reducing elastic energy throughout the whole solid. The effect of such a 
term on the transition from stable to unstable surface dynamics is depicted in the right- 
hand column of Fig. 4. As usual, we simplify the analysis by assuming that max(i?g) is 
achieved for q = {qmax, 0) and solve the two equations: (i) Rg = and (ii) dRg/dq = 0, from 
which we get S^f^{e) = S* = M^/AB.^^ > and \qmax\ = M/2B,^,^ > at the transition. 

In this analysis we have implicitly assumed that the transition from stable to unstable 
dynamics is "supercritical" - namely, that it is triggered by infinitesimal perturbations, and 
thus associated with a change of sign of max(i?g). It is also possible that the transition is 
"subcritical", and occurs at parameters for which the linear stability analysis, Eq. ( |3T1) yields 
max(i?q) < 0. If the transition is subcritical, then the characteristic wavelength may not 
diverge even if the linear dispersion is of the form (13T1) . It is possible to discern supercritical 
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from subcritical transitions by probing signatures of hysteretic behavior (associated with 
subcritical but not with supercritical transitions), and by carefully analyzing the kinetics of 
pattern formation. A necessary condition for the existence of a subcritical transition is that 
the leading nonlinear contributions to the dynamics have a destabilizing effect (unlike the 
stabilizing nonlinear terms derived in 2l|). Because our analysis is restricted to the linear 
dynamics we will not pursue this possibility further here. 



VI. CONCLUSIONS 



While the possibihty of producing patterned surfaces has attracted significant attention 
recently, few experiments have focused on regions in parameter space where dynamically 
stable, smooth surfaces are observed. The existence of these stable regions contradicts the 
Bradley-Harper stability analysis, but this is only part of the reason for their importance: we 
have argued in this paper that the emergence of stable surfaces provides important insights 
into the surface dynamics, that are critical for the development of a nonlinear theory of 
pattern formation in any parameters regime of ion sputtering. Our major messages are: 

1. The Bradley-Harper prediction regarding the instability of ion-bombarded surfaces 
to perpendicular mode ripples follows from a broad class of purely erosive response 
functions. This robustness may explain why the Bradley-Harper picture seems to 
describe correctly many observations of pattern evolution on ion sputtered surfaces. 

2. Various types of non-erosive response can change the sign of the coefficient of the 
second spatial derivative and thereby change the stability of surfaces to the emergence 
of large scale patterns. In particular, modifications of the response can lead to linear 
stability of smooth surfaces at various ranges of beam angles. These changes can 
be accompanied by no observable modification of the yield curve. Evidence for such 
modifications should thus come from atomistic simulations or from experiments that 
are capable of probing the local surface response to a single ion impact. 

3. Careful analysis of qualitative features of the pattern near the transition between sta- 
bility and instability of a fiat surface, in particular the existence or lack of divergence of 
the pattern wavelength at the transition, enable us to determine conclusively whether 
nonlocal mechanisms significantly affect the surface dynamics. The outcome of this 
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analysis is extremely important: because the existence of nonlocal terms qualitatively 
changes the linear dispersion relation, they must be included in the surface dynamics, 
even away from the transition regime. 

This paper focused on the linear dynamics of ion sputtered surfaces. In order to predict 
and control the fully developed patterns it is necessary to extend this to a nonlinear analysis. 
The existence of a stable-unstable transition at a critical beam angle 9c presents an excellent 
opportunity for quantitative predictions about pattern formation. Typically, near such a 
transition only a few Fourier modes are unstable, and the morphology of evolving patterns 
can generally be described by a weakly nonlinear "amplitude equation", whose form is 



universal and is determined almost solely by symmetry considerations [38|] . In other contexts, 
such amplitude equations have been enormously successful at predicting the shape of the 
selected patterns and many more features of their dynamics. Such an approach has not been 
tried so far for ion sputtered surfaces, apparently because it has been assumed that there is 
no continuous control parameter whose variation may change the stability of fiat surfaces. 
Recognizing that the beam angle is exactly such a parameter, at least for certain surfaces 
and ion types and energies, may enable the application of this invaluable theoretical tool to 
quantitative study of pattern formation on ion sputtered surfaces. 

We hope that the theoretical directions outlined in this paper will trigger experimental 
and computational work that will lead to better understanding of the surface response to 
ion impact and its relevance to large scale surface dynamics, and to better characterization 
of the transition from stability to instability of flat surfaces. We believe that such insights 
will be important to the development of a quantitative theory that will predict whether and 
what types of patterns are formed on a sputtered surface for a given set of material and ion 
beam parameters. 
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